#!/bin/csh -f
##############################################################################
#	Script for tAo graphic output with GMT 4.0
#	Daniel Garcia-Castellanos
##############################################################################
#Syntax:  tao.gmt.job <prjname> <xmin> <xmax> <ymin> <ymax> <switch_sea> 
#			<switch_rheo>
##############################################################################
#Note that:
# -This script uses bc unix command as a calculator. 
# -<xmin> <xmax> must be inside x range defined in .PRM parameters file.
##############################################################################

#Setting up variables:

set prj	= $argv[1]

source $tao_dir/script/tao.common.gmt.job


#1st PAGE of Postscript:

#
#MAIN GRAPHIC VERTICALLY EXAGGERATED:
set vertical_exag = 20
set graph_width = 15
#set ymincropkm = -17
#set ymincrop 	= `echo $ymincropkm \* 1000 | bc -l` 
set txt = $prj
set graph_height = `echo $vertical_exag \* $graph_width \* \($zmaxkm \- $zminkm\)\/\($xmax \- $xmin\) | bc -l` 
psbasemap -JX$graph_width/$graph_height -R$xmin/$xmax/$zminkm/$zmaxkm \
	-Ba50f10:"distance (km)":/2f.2:"elevation (km)":nSeW \
	-P -K -X3.5 -Y3 -U"tAo graphic output. Project: "$txt 		>! $ps
psbasemap -JX -R$xmin/$xmax/$zmin/$ymax -K -O -B::		>> $ps

#Draws Sea:
if ($switch_sea == 1) psxy $prj.sea.tao.gmt.tmp -JX -B/g1000 -R -O -K -L -G$color_sea	>> $ps

#Draws basement:
psxy $prj.basement.xz.tao.gmt.tmp -JX -R -O -K -L -G$color_heavy_basement 	>> $ps

#set uplift = `awk -v tn=$Timenow 'BEGIN{u=(tn+40)*3000/40; if (0) u=30*3000/40; print u; exit;}'`
#echo uplift = $uplift
#psxy -JX -R -W4ta -O -K << END >> $ps
#	-50	0
#	-50	$uplift
#	50	$uplift
#	50	0
#END

set col = 3
while ($col <= $numcols)
	awk '{print $1, $col}' col=$col $prj.pfl > $prj.0.tao.gmt.tmp
	awk '{print $1, $(col-1)}' col=$col $prj.2.tao.gmt.tmp >> $prj.0.tao.gmt.tmp

	if ($dens[$col] > 2845) then
		set color = $color_heavy_basement
	else
		if ($dens[$col] > 2550) then
			set color = $color_basement
		else
			if ($dens[$col] > 2400) then
				set color = $color_heavy_sediment
			else
				if ($dens[$col] > 1100) then
					set colant = `echo $col - 1 | bc`
					set age_cond = `awk -v a=$age[$colant] -v agelimit=$sedim_age_limit 'BEGIN{age=a*1.0; agelimit*=1.; if (age>=agelimit) print 1; else print 0.0}' `
					if ($age_cond) then
						set color = $color_sediment
					else
						set color = $color_old_sediment
					endif
				else
					if ($dens[$col] > 970) then
						set color = $color_sea
					else
						set color = $color_ice
					endif
				endif
			endif
		endif
	endif

	psxy $prj.0.tao.gmt.tmp -JX -R -O -K -L -M -W1 -G$color 	>> $ps
	set col = `echo $col + 1 | bc`
end

#if (-r $prj.CMP)  psxy $prj.CMP -JX -R -O -M -W10/255 -K		>> $ps
if (-r $prj.CMP)  psxy $prj.CMP -JX -R -O -M -W4ta/150/0/0 -K	>> $ps
if (-r $prj.FIT)  psxy $prj.FIT -JX -R -O -Ey -W1/255/150/150 -K  >> $ps
pstext	-JX -R -N -O -K -W255 << END >> $ps
	$xmin	$ymax	15 0 0 LT 	\ t=$Timenow My
END
psbasemap -JX -R -B:: -O -K			>> $ps


#
#TRUE SCALE GRAPHIC:
set vertshift = `echo $graph_height + 1.8 | bc -l`
set vertical_exag = 1
set graph_height = `echo $vertical_exag \* $graph_width \* \($ymaxkm \- $yminkm\)\/\($xmax \- $xmin\) | bc -l` 
psbasemap -JX$graph_width/$graph_height -R$xmin/$xmax/$yminkm/$ymaxkm \
	-Ba50f10/5f1:"elevation (km)":nSeW \
	-O -K -Y$vertshift 		>> $ps
psbasemap -JX -R$xmin/$xmax/$ymin/$ymax -K -O -B::		>> $ps

#Draws Sea:
if ($switch_sea == 1) psxy $prj.sea.tao.gmt.tmp -JX -B/g1000 -R -O -K -L -G$color_sea	>> $ps

#Draws basement:
psxy $prj.basement.xz.tao.gmt.tmp -JX -R -O -K -L -G$color_basement 	>> $ps

set col = 3
while ($col <= $numcols)
	awk '{print $1, $col}' col=$col $prj.pfl > $prj.0.tao.gmt.tmp
	awk '{print $1, $(col-1)}' col=$col $prj.2.tao.gmt.tmp >> $prj.0.tao.gmt.tmp

	if ($dens[$col] > 2895) then
		set color = $color_heavy_basement
	else
		if ($dens[$col] > 2550) then
			set color = $color_basement
		else
			if ($dens[$col] > 2400) then
				set color = $color_heavy_sediment
			else
				if ($dens[$col] > 1100) then
					set colant = `echo $col - 1 | bc`
					set age_cond = `awk -v a=$age[$colant] -v agelimit=$sedim_age_limit 'BEGIN{age=a*1.0; agelimit*=1.; if (age>agelimit) print 1; else print 0.0}' `
					if ($age_cond) then
						set color = $color_sediment
					else
						set color = $color_old_sediment
					endif
				else
					if ($dens[$col] > 970) then
						set color = $color_sea
					else
						set color = $color_ice
					endif
				endif
			endif
		endif
	endif

	psxy $prj.0.tao.gmt.tmp -JX -R -O -K -L -M -W1 -G$color 	>> $ps
	set col = `echo $col + 1 | bc`
end

#if (-r $prj.CMP)  psxy $prj.CMP -JX -R -O -M -W10/255 -K		>> $ps
if (-r $prj.CMP)  psxy $prj.CMP -JX -R -O -M -W4ta/150/0/0 -K	>> $ps
if (-r $prj.FIT)  psxy $prj.FIT -JX -R -O -Ey -W1/255/150/150 -K  >> $ps
psbasemap -JX -R -B:: -O -K			>> $ps



#
#PARAMETERS PROFILE:
set graph_height = 2.5
if (-r $prj.eros) then
set vertshift = `echo $graph_height + .9 | bc -l`

awk '{	if (NR==2) for (i=3;i<=NF;i++) {dens[i]=$i;}\
	if (NR>2){ thick=0;\
		   for (i=3;i<=NF;i++) if (dens[i]>2400) thick+=($i-$(i-1)); \
		   print $1, thick;\
		  }\
}' $prj.pfl > $prj.basam_thick.tao.gmt.tmp

awk '{	if (NR==2) for (i=3;i<=NF;i++) {dens[i]=$i;}\
	if (NR>2){ thick=0;\
		   for (i=3;i<=NF;i++) {if (dens[i]<=2400 && dens[i]>1100) thick+=($i-$(i-1)); } \
		   print $1, thick/1e3;\
		  }\
}' $prj.pfl > $prj.sed_thick.tao.gmt.tmp

awk '{print $1, $2/1e3}'     $prj.eros | psxy -JX$graph_width/$graph_height      -R$xmin/$xmax/-2/8 -W4/150/100/0 -Ba50f10/2f1g10:"total eros./sed. (km)":NseW -K -Y$vertshift -O >> $ps
#paste $prj.basam_thick_ref $prj.basam_thick.tao.gmt.tmp | awk '{print $1,$4-$2}' | \
#	psxy -JX -R -O -K -W3/255/0/0 >> $ps 
#psxy $prj.sed_thick.tao.gmt.tmp -JX -R -O -K -W3/200/170/0 >> $ps 

set vertshift = `echo $graph_height + .9 | bc -l`
gmtset BASEMAP_FRAME_RGB +0/0/0
psbasemap -JX$graph_width/$graph_height -R$xmin/$xmax/-10/10 -Ba50f10:"":/Ns -Y$vertshift -O -K  >> $ps
gmtset BASEMAP_FRAME_RGB +180/0/0;     awk '{print $1, $3/1e3}' $prj.eros | psxy -JX -R$xmin/$xmax/-1/1 -W4/180/0/0 -O -K -B/a1f1g100:"erosion rate (mm/yr)":W >> $ps
gmtset BASEMAP_FRAME_RGB +150/100/0;   awk '{print $1, $6}'     $prj.eros | psxy -JX -R$xmin/$xmax/0/.002 -W4/150/100/0 -O -K -B/.001:"sed. load (kg/s/m)":E	>> $ps
gmtset BASEMAP_FRAME_RGB +0/0/0; 
set total_rain = `awk '(NR>3){total+=$8}END{printf("%.0f", total)}' $prj.eros`
if ($total_rain) then 
psbasemap -JX$graph_width/$graph_height -R$xmin/$xmax/-10/10 -Ba50f10:"":/Ns -Y$vertshift -O -K >> $ps
gmtset BASEMAP_FRAME_RGB +0/0/255;     awk '{print $1, $5}'     $prj.eros | psxy -JX -R$xmin/$xmax/0/.003 -W4/0/0/255 -O -K -B/.001:"water disch. (m2/s)":E	>> $ps
gmtset BASEMAP_FRAME_RGB +0/150/255;   awk '{print $1, $8/1e3}'     $prj.eros | psxy -JX -R$xmin/$xmax/0/4 -W4/0/150/255 -O -K -B/a2f1:"P, E (m/yr)":W	>> $ps
gmtset BASEMAP_FRAME_RGB +100/180/255; awk '{print $1, $9/1e3}'     $prj.eros | psxy -JX -R$xmin/$xmax/0/4 -W6/100/180/255ta -O -K	>> $ps
gmtset BASEMAP_FRAME_RGB +0/0/0; 
endif
endif

if (-r $prj.xzt) then
gmtset BASEMAP_FRAME_RGB +0/0/0;   
psbasemap -JX$graph_width/$graph_height -R$xmin/$xmax/-1.5/.5 -Ba50f10:"":/Ns -Y$vertshift -O -K >> $ps
gmtset BASEMAP_FRAME_RGB +0/0/0;     awk '{print $1, -$(NF-1)/1e3}'     $prj.xzt | psxy -JX -R -W4/0/0/0 -O -K -B/.5:"vert. deflection (km)":W	>> $ps
endif

psbasemap -JX -R -O -Bnsew			>> $ps








#2nd PAGE of Postscript:
#Gravity Anomaly graphic:
if (-r $prj.xg) then 
	awk '{if (NR>=3) print $1, $2}' $prj.xg \
		| psxy -JX$graph_width/10 -R$xmin/$xmax/-200/200 -W4/155/0/0 -Y2.7 -X3 -K \
		-B:.$prj\:a50f5g25:"distance (km)":/100f10g50:"Gravity Anomaly (mGal)":nSeW\
		>> $ps
	awk '{if (NR>=3) print $1, $3}' $prj.xg \
		| psxy -JX -R$xmin/$xmax/-20/20 -W4/0/155/0 -O \
		-B/10f1:"Geoid Height (m)":E \
		>> $ps
	if (-r $prj.xg.CMP) then 
		awk '{print $1, $2}' $prj.xg.CMP \
			| psxy -JX -R -W4/155/0/0 -St.1 -O -K \
			>> $ps
		awk '{print $1, $3}' $prj.xg.CMP \
			| psxy -JX -R -W4/155/0/0 -St.1 -O -K \
			>> $ps
	endif
endif


if ($switch_rheo == 0) then 
	rm -f  $prj.*.tao.gmt.tmp
	exit
endif

#3rd PAGE of Postscript:
#
#DEFLECTION & STRESSES
if (-r $prj.strs) then
	#rm -f $prj.temp.grd $prj.strs.grd 
	#awk '{print $1,$2,$4}' $prj.temp | surface -G$prj.temp.grd -I2/1 -R$xmin/$xmax/-5/70 -H2 
	surface	$prj.strs -G$prj.strs.grd -I2/1 -R$xmin/$xmax/-2/70 -H2 
	#if (-r $prj.temp.grd) rm -f $prj.temp
	#if (-r $prj.strs.grd) rm -f $prj.strs
endif
psbasemap -JX$graph_width/-5 -R$xmin/$xmax/-2/70 -Ba50f10:"distance (km)":/20f10:"Depth (km)":nSeW \
	-X4 -Y3 -K 										>> $ps
psclip	$prj.stress.clip.tao.gmt.tmp -JX -R -O -K 						>> $ps
grdimage	$prj.strs.grd -C$prj.stress_cpt.tao.gmt.tmp -JX -R -O -K  			>> $ps
grdcontour	$prj.strs.grd -JX -R -C$prj.contours1.z.tao.gmt.tmp -A10f10 -W1/255 -K -O 	>> $ps
#grdcontour	$prj.strs.grd -JX -R -C$prj.contours2.z.tao.gmt.tmp -A10f10 -W1/255 -K -O 	>> $ps
psclip -C -O -K 										>> $ps
psxy	$prj.stress.clip.tao.gmt.tmp -JX -R -O -K -W7						>> $ps
psbasemap -JX -R -Bnsew -O -K 	>> $ps

psscale -C$prj.stress_cpt.tao.gmt.tmp -B:."stress (MPa)": -D5/-.7/7/.2h -O -K			>> $ps

if (-r $prj.xzt) awk '{print($1,-$(NF-1))}' $prj.xzt | \
	psxy -JX$graph_width/3.75 -R$region -Ba50f10/2000:"deflection (m)":neW \
	-W7 -H2 -Y5.5 -O -K  	>> $ps


psbasemap -JX$graph_width/-2 -R$xmin/$xmax/0/70 -B:.$prj\:a50f10/20f10:"thickness (km)":NeW \
	-Y4.5 -O -K	>> $ps
if (-r $prj.CRUST)  awk '{print $1/1000, $2}' $prj.CRUST |psxy  -JX -R$xmin/$xmax/0/60000 -W10/60 -L -G170 -O -K  	>> $ps
if (-r $prj.UCRUST) awk '{print $1/1000, $2}' $prj.UCRUST |psxy  -JX -R -O -K -L -G200 -W10/130 	>> $ps
psxy $prj.eeth -JX -R$xmin/$xmax/0/60000 -O -W10 -K -H1 					>> $ps
pstext	-JX	-R  -O 									<< END	>> $ps
#	0	25000	13  0   1  1	Te
#	-30	48000	11  0   1  2	H
END



#4th PAGE of Postscript:
#MAXIMUM MOMENT POINT:
if (-r $prj.ysen) then
psxy $prj.ysen -JX12.5/-22 -R-1e3/1e3/0/1e2 -W8 -P -K -: \
	-B:."Yield Stress Envelope '$prj'":a200f100:"Stress (MPa)":/a10f5:"Depth (km)":NsEW 	>> $ps
awk '{print($1,$3)}' $prj.ysen | psxy -JX -R -W8 -O -K -:  					>> $ps
awk '{print($1,$5)}' $prj.ysen | psxy -JX -R -W5/140 -H1 -M -O -K -:				>> $ps
awk '{print($1,$4)}' $prj.ysen | psxy -JX -R -W12to -B200f100:"Temperature (C)":nSew -O -:  	>> $ps
endif

rm -f  $prj.*.tao.gmt.tmp
